extern "C" {
#include "mcce.h"
}

#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <set>
#include <string>
using namespace std;

const string PDBF = "step2_out.pdb";
const string HBOUT = "hb.dat";
const float DFAR = 3.2;       // the distance used to determine a hydrogen bond
const float DNEAR = 1.2;       // 1.2 < d < 3.2

int is_hb(CONF *conf1, CONF *conf2);
int load_head3lst();

FILE *h_fp;
static RES conflist;
const float PI = 3.1415926;

int main()
{
	db_open();
	if (init()) {
		db_close();
		printf("Help message: double check file \"run.prm\" in current directory.\n");
		return USERERR;
	}

	FILE *fp;
	PROT prot;
	if (!(fp=fopen(STEP2_OUT, "r"))) {
		printf("   \"No step 2 output \"%s\".\n", STEP2_OUT);
		return USERERR;
	}
	prot = load_pdb(fp);
	if (prot.n_res == 0) {
		printf("   There are errors in pdb file, quiting ...\n");
		return USERERR;
	}

	id_conf(prot);
	get_connect12(prot);

	load_head3lst();
	int n_conf = conflist.n_conf;

	int i, j, k;
	//    for (k=0; k<n_conf; k++) printf("%s\t%d\n", conflist.conf[k].uniqID, conflist.conf[k].iConf);

	for (i=0; i<prot.n_res; i++) {
		for (j=1; j<prot.res[i].n_conf; j++) {
			strncpy(prot.res[i].conf[j].confName, prot.res[i].conf[j].uniqID, 5);
			prot.res[i].conf[j].confName[5] = '\0';
			// adjust the conf id to avoid id diff between step2 and head3.lst
			for (k=0; k<n_conf; k++) {
				if (!strcmp(conflist.conf[k].uniqID, prot.res[i].conf[j].uniqID)) {
					prot.res[i].conf[j].iConf = conflist.conf[k].iConf;
					break;
				}
			}
		}
	}



	h_fp = fopen("hah.txt", "w");

	// hb_fp is a binary file to store the hbpw matrix, the first 4 bytes give the conformer number n_conf
	FILE *hb_fp = fopen(HBOUT.c_str(), "wb");
	fwrite(&n_conf, sizeof(int), 1, hb_fp);

	int i_res, j_res, i_conf, j_conf;

	// hbpw is a matrix to store the hb connection between each two conf, the elem is 0 or 1
	FILE *fp_res = fopen("reshbond.txt", "w");
	char *resHb = (char *) calloc(prot.n_res * prot.n_res, sizeof(char));

	char *hbpw = (char *) calloc(n_conf * n_conf, sizeof(char));
	for (i_res=0; i_res<prot.n_res; i_res++) {
		for (i_conf=1; i_conf<prot.res[i_res].n_conf; i_conf++) {
			for (j_res=0; j_res<prot.n_res; j_res++) {
				if (j_res == i_res) continue;
				for (j_conf=1; j_conf<prot.res[j_res].n_conf; j_conf++) {

					/* there is hbond from donor i_conf to accepter j_conf */
					if (is_hb(&prot.res[i_res].conf[i_conf], &prot.res[j_res].conf[j_conf])) {
						hbpw[prot.res[i_res].conf[i_conf].iConf * n_conf + prot.res[j_res].conf[j_conf].iConf] = 1;
						resHb[i_res * prot.n_res + j_res] = 1;
					}
				}
			}
		}
	}
	/*
    for (i=0; i<n_conf; i++) {
        for (j=0; j<i; j++) {
            hbpw[i * n_conf + j] = hbpw[j * n_conf + i];
        }
    }
	 */
	fwrite(hbpw, sizeof(char), n_conf * n_conf, hb_fp);
	fclose(hb_fp);

	set <int> resInHbNet;
	for (i_res=0; i_res<prot.n_res; i_res++) {
		for (j_res=0; j_res<prot.n_res; j_res++) {
			if (j_res == i_res) continue;
			if (resHb[i_res * prot.n_res + j_res] == 1) {
				fprintf(fp_res, "%s%c%04d\t%s%c%04d\n", prot.res[i_res].resName, prot.res[i_res].chainID, prot.res[i_res].resSeq,
						prot.res[j_res].resName, prot.res[j_res].chainID, prot.res[j_res].resSeq);
				resInHbNet.insert(i_res);
				resInHbNet.insert(j_res);
			}
		}
	}
	fclose(fp_res);

	FILE *fp_resInHbNet = fopen("resInHbNet.txt", "w");
	for (i_res=0; i_res<prot.n_res; i_res++) {
		if (resInHbNet.find(i_res) != resInHbNet.end()) {
			fprintf(fp_resInHbNet, "%s%c%04d\n", prot.res[i_res].resName, prot.res[i_res].chainID, prot.res[i_res].resSeq);
		}
	}
	fclose(fp_resInHbNet);

	printf("n_res: %d\n", prot.n_res);
	db_close();
	return 0;
}

int is_hb(CONF *conf1, CONF *conf2)
{
	STRINGS Datoms, Aatoms;
	int iD, iA, Dseq, Aseq;
	float d=0.0;
	float ang = 0.0;
	int isHb = 0;
	/** To establish hydrogen bond between two conformers, the following criteria must be met:
	 * 1. one conformer is h donor, the other is acceptor.
	 * 2. the distance between one donor atom and one acceptor atom is between DNEAR and DFAR.
	 * 3. the angle of h bond should be no less than 90.
	 */
	if (!param_get("HDONOR", conf1->confName, "", &Datoms)) {
		if (!param_get("HACCEPT", conf2->confName, "", &Aatoms)) {
			for (iD=0; iD<Datoms.n; iD++) {
				param_get("IATOM", conf1->confName, Datoms.strings[iD], &Dseq);
//				if (!strcmp(conf1->confName, "HOH-1")) {
//					printf("Donor: %s, Datoms.n: %d, Dseq: %d, atom name %s\n", conf1->confName,
//						Datoms.n, Dseq, Datoms.strings[iD]);
//				}
				for (iA=0; iA<Aatoms.n; iA++) {
					param_get("IATOM", conf2->confName, Aatoms.strings[iA], &Aseq);
					d = ddvv(conf1->atom[Dseq].xyz, conf2->atom[Aseq].xyz);

//					d = ddvv(conf1->atom[Dseq].connect12[0]->xyz, conf2->atom[Aseq].xyz);
					if (d > DNEAR * DNEAR && d < DFAR * DFAR) {
						ang = avv(vector_vminusv((conf1->atom[Dseq].connect12[0])->xyz, conf1->atom[Dseq].xyz),
								vector_vminusv(conf2->atom[Aseq].xyz, conf1->atom[Dseq].xyz)) * 180.0 / PI;
						if (abs(ang) < 90.0) continue;
						fprintf(h_fp, "%s\t%s\t%s~%s--%s\t%.2f\t%.0f\n", conf1->uniqID, conf2->uniqID,
								(conf1->atom[Dseq].connect12[0])->name, conf1->atom[Dseq].name, conf2->atom[Aseq].name, sqrt(d), ang);
						isHb = 1;
					}
				}
			}
		}
	}

	return isHb;
}

int load_head3lst()
{
	FILE *fp;
	char sbuff[MAXCHAR_LINE];
	char stemp[MAXCHAR_LINE];
	CONF conf_temp;
	char notfound;
	int iconf, ires;
	int kr;
	int counter;

	conflist.n_conf = 0;
	conflist.conf   = NULL;

	if (!(fp=fopen(FN_CONFLIST3, "r"))) {
		printf("   FATAL: Can't open file %s\n", FN_CONFLIST3);
		return USERERR;
	}
	fgets(sbuff, sizeof(sbuff), fp); /* skip the first line */
	counter = 0;
	while(fgets(sbuff, sizeof(sbuff), fp)) {
		/* load this line to a conf template */
		if (strlen(sbuff) < 20) continue;
		sscanf(sbuff, "%d %s %c %f %f %f %f %d %d %f %f %f %f %f %f %s", &conf_temp.iConf,
				conf_temp.uniqID,
				&conf_temp.on,
				&conf_temp.occ,
				&conf_temp.netcrg,
				&conf_temp.Em,
				&conf_temp.pKa,
				&conf_temp.e,
				&conf_temp.H,
				&conf_temp.E_vdw0,
				&conf_temp.E_vdw1,
				&conf_temp.E_tors,
				&conf_temp.E_epol,
				&conf_temp.E_dsolv,
				&conf_temp.E_extra,
				conf_temp.history);

		conf_temp.E_TS = 0.0; /* initialize entropy effect at the time of loading conflist */

		/* rescale */
		conf_temp.E_vdw0 *= env.scale_vdw0;
		conf_temp.E_vdw1 *= env.scale_vdw1;
		conf_temp.E_epol *= env.scale_ele;
		conf_temp.E_tors *= env.scale_tor;
		conf_temp.E_dsolv*= env.scale_dsolv;

		strncpy(conf_temp.resName, conf_temp.uniqID, 3); conf_temp.resName[3] = '\0';
		strncpy(conf_temp.confName, conf_temp.uniqID, 5); conf_temp.confName[5] = '\0';
		conf_temp.chainID = conf_temp.uniqID[5];
		strncpy(stemp, conf_temp.uniqID+6, 4); stemp[4] = '\0';
		conf_temp.resSeq = atoi(stemp);
		conf_temp.iCode = conf_temp.uniqID[10];
		conf_temp.n_atom = 0;
		if (conf_temp.on == 't' || conf_temp.on == 'T') conf_temp.on = 't';
		else conf_temp.on = 'f';
		conf_temp.iConf = counter;
		/* creating conflist */
		iconf = ins_conf(&conflist, conflist.n_conf, 0);
		cpy_conf(&conflist.conf[iconf], &conf_temp);
		counter++;
	}

	return 0;
}
